#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <assert.h>
#include <Python.h>
#include <numpy/arrayobject.h>

//header structure
struct io_header
{
  int      npart[6];
  double   mass[6];
  double   time;
  double   redshift;
  int      flag_sfr;
  int      flag_feedback;
  int      npartTotal[6];
  int      flag_cooling;
  int      num_files;
  double   BoxSize;
  double   Omega0;
  double   OmegaLambda;
  double   HubbleParam; 
  int      FlagAge;
  int      FlagMetals;
  int      NallHW[6];
  int      flag_entr_ics;
  char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8 - 2*4 - 6*4 - 4];  /* fills to 256 Bytes */
};

//read positions
static PyObject *
readone(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,pos,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  
  //ask if args is a string and get position wanted (x,y,z,vx,vy,vz) = (0,1,2,3,4,5)
  if (!PyArg_ParseTuple(args, "sii", &file, &nfiles, &pos))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readOne syntax is gadgetPyIO.readOne(<filename>,<nfiles>,pos) where pos is one of (0,1,2,3,4,5) = (x,y,z,vx,vy,vz)");
      return NULL;
    }
  if(pos < 0 || pos > 5) //error checking pos
    {
      PyErr_Format(PyExc_IndexError,"gadgetPyIO.readOne needs the phase space coordinate to match one of (x,y,z,vx,vy,vz) = (0,1,2,3,4,5)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_FLOAT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      if(pos < 3) //read a position
	{
	  fread(&dummy, sizeof(dummy), 1, fd);
	  for(k=0,pc_new=pc;k<6;k++)
	    {
	      for(n=0;n<header.npart[k];n++)
		{
		  switch(pos)
		    {
		    case 0:
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      fseek(fd,2*sizeof(float),SEEK_CUR);
		      break;
		    case 1:
		      fseek(fd,sizeof(float),SEEK_CUR);
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      fseek(fd,sizeof(float),SEEK_CUR);
		      break;
		    case 2:
		      fseek(fd,2*sizeof(float),SEEK_CUR);
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      break;
		    }
		  pc_new++;
		}
	    }
	  fread(&dummy, sizeof(dummy), 1, fd);
	}
      else  //read a velocity
	{
	  //skip positions
	  fread(&dummy, sizeof(dummy), 1, fd);
	  fseek(fd,dummy,SEEK_CUR);
	  fread(&dummy, sizeof(dummy), 1, fd);
	  
	  fread(&dummy, sizeof(dummy), 1, fd);
	  for(k=0,pc_new=pc;k<6;k++)
	    {
	      for(n=0;n<header.npart[k];n++)
		{
		  switch(pos)
		    {
		    case 3:
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      fseek(fd,2*sizeof(float),SEEK_CUR);
		      break;
		    case 4:
		      fseek(fd,sizeof(float),SEEK_CUR);
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      fseek(fd,sizeof(float),SEEK_CUR);
		      break;
		    case 5:
		      fseek(fd,2*sizeof(float),SEEK_CUR);
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      break;
		    }
		  pc_new++;
		}
	    }
	  fread(&dummy, sizeof(dummy), 1, fd);
	}
      
      fclose(fd);
    }
  
  array = (PyArrayObject *) PyArray_Return(array);
  
  return (PyObject *) array;
}

//read ids
static PyObject *
readids(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readIds syntax is gadgetPyIO.readIds(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_INT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      //skip positions
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //skip velocities
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);

      //get ids!
      fread(&dummy, sizeof(dummy), 1, fd);
      for(k=0,pc_new=pc;k<6;k++)
	{
	  for(n=0;n<header.npart[k];n++)
	    {
	      fread(PyArray_GETPTR1(array,pc_new),sizeof(int),1,fd);
	      pc_new++;
	    }
	}
      fread(&dummy, sizeof(dummy), 1, fd);
      
      fclose(fd);
    }
  
  array = (PyArrayObject *) PyArray_Return(array);
  
  return (PyObject *) array;
}

//read Masses
static PyObject *
readmasses(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readMasses syntax is gadgetPyIO.readMasses(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_FLOAT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      //skip positions
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //skip velocities
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);

      //skip ids
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //read masses
      if(ntot_withmasses>0)
	fread(&dummy, sizeof(dummy), 1, fd);
      for(k=0,pc_new=pc;k<6;k++)
        {
          for(n=0;n<header.npart[k];n++)
            {
              if(header.mass[k] == 0)
		fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
	      else
		*((float*)PyArray_GETPTR1(array,pc_new)) = header.mass[k];
	      pc_new++;
	    }
        }
      
      fclose(fd);
    }
  
  array = (PyArrayObject *) PyArray_Return(array);
  
  return (PyObject *) array;
}


//read header
static PyObject *
printHeaders(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.printHeaders syntax is gadgetPyIO.printHeaders(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
            
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  fprintf(stdout,"total number of particles = %d\n",NumPart);
	  fprintf(stdout,"total number of gas particles = %d\n",Ngas);
	}
      fprintf(stdout,"Header #%d:\n",i);
      fprintf(stdout,"\tnpart[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.npart[0],header.npart[1],header.npart[2],header.npart[3],header.npart[4],header.npart[5]);
      fprintf(stdout,"\tmass[0-5] = (%e|%e|%e|%e|%e|%e)\n",header.mass[0],header.mass[1],header.mass[2],header.mass[3],header.mass[4],header.mass[5]);
      fprintf(stdout,"\ttime = %e\n",header.time);
      fprintf(stdout,"\tredshift = %e\n",header.redshift);
      fprintf(stdout,"\tflag_sfr = %d\n",header.flag_sfr);
      fprintf(stdout,"\tflag_feedback = %d\n",header.flag_feedback);
      fprintf(stdout,"\tnpartTotal[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.npartTotal[0],header.npartTotal[1],header.npartTotal[2],
	      header.npartTotal[3],header.npartTotal[4],header.npartTotal[5]);
      fprintf(stdout,"\tflag_cooling = %d\n",header.flag_cooling);
      fprintf(stdout,"\tnum_files = %d\n",header.num_files);
      fprintf(stdout,"\tBoxSize = %e\n",header.BoxSize);
      fprintf(stdout,"\tOmega0 = %f\n",header.Omega0);
      fprintf(stdout,"\tOmegaLambda = %f\n",header.OmegaLambda);
      fprintf(stdout,"\tHubbleParam = %f\n",header.HubbleParam);
      fprintf(stdout,"\tFlagAge = %d\n",header.FlagAge);
      fprintf(stdout,"\tFlagMetals = %d\n",header.FlagMetals);
      fprintf(stdout,"\tNallHW[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.NallHW[0],header.NallHW[1],header.NallHW[2],header.NallHW[3],header.NallHW[4],header.NallHW[5]);
      fprintf(stdout,"\tflag_entr_ics = %d\n",header.flag_entr_ics);
      
      fclose(fd);
    }
  
  Py_RETURN_NONE;
}

//read header
static PyObject *
printHeadersTypeTwo(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  char fileblockname[5]={"    "};
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.printHeaders syntax is gadgetPyIO.printHeaders(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy,sizeof(dummy),1,fd);
      fread(&fileblockname,4*sizeof(char),1,fd);
      fseek(fd,dummy-4*sizeof(char),SEEK_CUR);
      fread(&dummy,sizeof(dummy),1,fd);
      
      fread(&dummy,sizeof(dummy),1,fd);
      fread(&header,sizeof(header),1,fd);
      fread(&dummy,sizeof(dummy),1,fd);
            
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  fprintf(stdout,"total number of particles = %d\n",NumPart);
	  fprintf(stdout,"total number of gas particles = %d\n",Ngas);
	}
      fprintf(stdout,"Header #%d:\n",i);
      fprintf(stdout,"\tnpart[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.npart[0],header.npart[1],header.npart[2],header.npart[3],header.npart[4],header.npart[5]);
      fprintf(stdout,"\tmass[0-5] = (%e|%e|%e|%e|%e|%e)\n",header.mass[0],header.mass[1],header.mass[2],header.mass[3],header.mass[4],header.mass[5]);
      fprintf(stdout,"\ttime = %e\n",header.time);
      fprintf(stdout,"\tredshift = %e\n",header.redshift);
      fprintf(stdout,"\tflag_sfr = %d\n",header.flag_sfr);
      fprintf(stdout,"\tflag_feedback = %d\n",header.flag_feedback);
      fprintf(stdout,"\tnpartTotal[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.npartTotal[0],header.npartTotal[1],header.npartTotal[2],
	      header.npartTotal[3],header.npartTotal[4],header.npartTotal[5]);
      fprintf(stdout,"\tflag_cooling = %d\n",header.flag_cooling);
      fprintf(stdout,"\tnum_files = %d\n",header.num_files);
      fprintf(stdout,"\tBoxSize = %e\n",header.BoxSize);
      fprintf(stdout,"\tOmega0 = %f\n",header.Omega0);
      fprintf(stdout,"\tOmegaLambda = %f\n",header.OmegaLambda);
      fprintf(stdout,"\tHubbleParam = %f\n",header.HubbleParam);
      fprintf(stdout,"\tFlagAge = %d\n",header.FlagAge);
      fprintf(stdout,"\tFlagMetals = %d\n",header.FlagMetals);
      fprintf(stdout,"\tNallHW[0-5] = (%d|%d|%d|%d|%d|%d)\n",header.NallHW[0],header.NallHW[1],header.NallHW[2],header.NallHW[3],header.NallHW[4],header.NallHW[5]);
      fprintf(stdout,"\tflag_entr_ics = %d\n",header.flag_entr_ics);
      
      fclose(fd);
    }
  
  Py_RETURN_NONE;
}

//read Types
static PyObject *
readparttype(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readTypes syntax is gadgetPyIO.readTypes(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_INT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      for(k=0,pc_new=pc;k<6;k++)
        {
          for(n=0;n<header.npart[k];n++)
            {
              *((int*)PyArray_GETPTR1(array,pc_new)) = k;
	      ++pc_new;
	    }
	}
      
      fclose(fd);
    }
  
  array = (PyArrayObject *) PyArray_Return(array);
  
  return (PyObject *) array;
}

//read Types
static PyObject *
readparttypetypetwo(PyObject *self, PyObject *args)
{
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses;
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  float *temp,tempfloat;
  struct io_header header;
  FILE *fd;
  char buf[500];
  char fileblockname[5];
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "si", &file, &nfiles))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readTypesTypeTwo syntax is gadgetPyIO.readTypesTypeTwo(<filename>,<nfiles>)");
      return NULL;
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy,sizeof(dummy),1,fd);
      fread(&fileblockname,4*sizeof(char),1,fd);
      fseek(fd,dummy-4*sizeof(char),SEEK_CUR);
      fread(&dummy,sizeof(dummy),1,fd);
      if(strcmp(fileblockname,"HEAD") == 0)
	{
	  fread(&dummy, sizeof(dummy), 1, fd);
	  fread(&header, sizeof(header), 1, fd);
	  fread(&dummy, sizeof(dummy), 1, fd);
	}
      else
	{
	  PyErr_Format(PyExc_IOError,"first block is not labeled 'HEAD'!");
	  fclose(fd);
	  return NULL;
	}
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_INT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      for(k=0,pc_new=pc;k<6;k++)
        {
          for(n=0;n<header.npart[k];n++)
            {
              *((int*)PyArray_GETPTR1(array,pc_new)) = k;
	      ++pc_new;
	    }
	}
      
      fclose(fd);
    }
  
  array = (PyArrayObject *) PyArray_Return(array);
  
  return (PyObject *) array;
}

//read gas properties
static PyObject *
readgasprops(PyObject *self, PyObject *args)
{
  //FLAG_NOGAS is used to indicate that we have found gas particles in the file
  //   - we will assume that we have no particles at first
  //   - if we find them, the we set FLAG_NOGAS = 0 and return the array
  //   - else we decrement the refernce to the array we made and return nothing
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses,FLAG_NOGAS=1; 
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph,DOCONV=0;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  double tempdouble;
  struct io_header header;
  FILE *fd;
  char buf[500];
  const char *blocknamein;
  char blockname[5]={"    "};
  
  double GRAVITY, BOLTZMANN, PROTONMASS;
  double UnitLength_in_cm, UnitMass_in_g, UnitVelocity_in_cm_per_s;
  double UnitTime_in_s, UnitDensity_in_cgs, UnitPressure_in_cgs, UnitEnergy_in_cgs;  
  double G, Xh;
  double MeanWeight, u, gamma;
 
  //physical constants in cgs units
  GRAVITY   = 6.672e-8;
  BOLTZMANN = 1.3806e-16;
  PROTONMASS = 1.6726e-24;
  
  //internal unit system of the code - "preferred" units in Gadget
  UnitLength_in_cm = 3.085678e21;   //code length unit in cm/h
  UnitMass_in_g = 1.989e43;         //code mass unit in g/h
  UnitVelocity_in_cm_per_s = 1.0e5; //code unit in cm/s
      
  //ask if args is a string
  if(!PyArg_ParseTuple(args, "sis|iddd", &file, &nfiles, &blocknamein, &DOCONV, &UnitLength_in_cm, &UnitMass_in_g, &UnitVelocity_in_cm_per_s))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readGasProps syntax is gadgetPyIO.readGasProps(<filename>,<nfiles>,'propstring',[DOCONV,UnitLength,UnitMass,UnitVelocity])\nSee the module doc string for the 'propstring' conventions.  (In ipython 'gadgetPyIO.readGasProps?' will print the doc string.)");
      return NULL;
    }
  
  //error check blockname to be the correct size
  if(strlen(blocknamein) > 4)
    {
      PyErr_Format(PyExc_TypeError,"'blockname' must be a four or less character string!");
      return NULL;
    }
  else
    {
      //put into blockname and make sure it is upper case
      for(i=0;i<strlen(blocknamein);++i)
	blockname[i] = blocknamein[i];
      for(i=strlen(blocknamein);i<4;++i)
	blockname[i] = ' ';
      for(i=0;i<4;++i)
	blockname[i] = toupper(blockname[i]);
    }
  
  //get the other units
  UnitTime_in_s = UnitLength_in_cm/UnitVelocity_in_cm_per_s;
  UnitDensity_in_cgs = UnitMass_in_g/pow(UnitLength_in_cm,3);
  UnitPressure_in_cgs = UnitMass_in_g/UnitLength_in_cm/pow(UnitTime_in_s,2);
  UnitEnergy_in_cgs = UnitMass_in_g*pow(UnitLength_in_cm,2)/pow(UnitTime_in_s,2);
  G = GRAVITY/pow(UnitLength_in_cm,3)*UnitMass_in_g*pow(UnitTime_in_s,2);
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s'\n",buf); 
      fflush(stdout);
      
      fread(&dummy, sizeof(dummy), 1, fd);
      fread(&header, sizeof(header), 1, fd);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_FLOAT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      //skip positions
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //skip velocities
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);

      //skip ids
      fread(&dummy, sizeof(dummy), 1, fd);
      fseek(fd,dummy,SEEK_CUR);
      fread(&dummy, sizeof(dummy), 1, fd);
      
      //skip masses if they exist
      if(ntot_withmasses>0)
        {
	  fread(&dummy, sizeof(dummy), 1, fd);
	  fseek(fd,dummy,SEEK_CUR);
	  fread(&dummy, sizeof(dummy), 1, fd);
	}
      
      //get gas props
      //they are stored in order
      //1) internal energy
      //2) density
      if(header.npart[0] > 0)
        {
	  FLAG_NOGAS = 0;
	  if(strcmp(blockname,"U   ") == 0)
	    {
	      //read energy
	      fread(&dummy, sizeof(dummy), 1, fd);
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  if(DOCONV)
		    {
		      tempdouble = ((double) (*((float*)PyArray_GETPTR1(array,pc_sph))))*UnitEnergy_in_cgs/UnitMass_in_g;
		      (*((float*)PyArray_GETPTR1(array,pc_sph))) = (float) tempdouble; 
		    }
		  pc_sph++;
		}
	      fread(&dummy, sizeof(dummy), 1, fd);
	    }
	  else if(strcmp(blockname,"RHO ") == 0)
	    {
	      //skip energy
	      fread(&dummy, sizeof(dummy), 1, fd);
	      fseek(fd,dummy,SEEK_CUR);
	      fread(&dummy, sizeof(dummy), 1, fd);
	      
	      //read density
	      fread(&dummy, sizeof(dummy), 1, fd);
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  if(DOCONV)
		    {
		      tempdouble = ((double) (*((float*)PyArray_GETPTR1(array,pc_sph))))*UnitDensity_in_cgs;
		      (*((float*)PyArray_GETPTR1(array,pc_sph))) = (float) tempdouble;
		    }
		  pc_sph++;
		}
	      fread(&dummy, sizeof(dummy), 1, fd);
	    }
	  else if(strcmp(blockname,"HSML") == 0)
	    {
	      //skip energy
	      fread(&dummy, sizeof(dummy), 1, fd);
	      fseek(fd,dummy,SEEK_CUR);
	      fread(&dummy, sizeof(dummy), 1, fd);
	      
	      //skip density
	      fread(&dummy, sizeof(dummy), 1, fd);
	      fseek(fd,dummy,SEEK_CUR);
	      fread(&dummy, sizeof(dummy), 1, fd);
	      
	      //read smoothing length
	      fread(&dummy, sizeof(dummy), 1, fd);
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  pc_sph++;
		}
	      fread(&dummy, sizeof(dummy), 1, fd);
	    }
	}
      else
	{
	  fprintf(stdout,"gadgetPyIO module found no gas particles in file %d!\n",i);
	}
      
      //advance pc_new for indexing
      for(k=0,pc_new=pc;k<6;k++)
	pc_new += header.npart[k];
      
      fclose(fd);
    }
  
  if(FLAG_NOGAS)
    {
      //not returning any gas info - decrement reference to avoid a memory leak!
      Py_DECREF(array);
      Py_RETURN_NONE;
    }
  else
    {
      array = (PyArrayObject *) PyArray_Return(array);
      return (PyObject *) array;
    }
}

//read gas properties
static PyObject *
readBlockTypeTwo(PyObject *self, PyObject *args)
{
  //FLAG_NOBLOCK is used to indicate that we have found the requested block in the snapshot
  //   - we will assume that we have not found it at first
  //   - if we find them, the we set FLAG_NOBLOCK = 0 and return the array
  //   - else we decrement the refernce to the array we made and return nothing
  const char *file;
  int i,j,k,nfiles,tempint,dummy,ntot_withmasses,FLAG_NOBLOCK=1; 
  int NumPart,Ngas,t,n,off,pc,pc_new,pc_sph;
  npy_intp dim[1] = {1};
  npy_intp ndim = 1;
  PyArrayObject *array;
  double tempdouble;
  struct io_header header;
  FILE *fd;
  char buf[500];
  const char *blocknamein;
  char blockname[5]={"    "};
  char fileblockname[5]={"    "};
  int pos=-1;
  
  //ask if args is a string
  if (!PyArg_ParseTuple(args, "sis|i", &file, &nfiles, &blocknamein, &pos))  
    {
      PyErr_Format(PyExc_TypeError,
		   "gadgetPyIO.readBlockTypeTwo syntax is gadgetPyIO.readGasProps(<filename>,<nfiles>,'blockname',[pos])\n'blockname' is the name of the Gadget2 block\nIf you specify a position, velocity, or acceleration block, you must specify which coordinate to read (x,y,z) = (0,1,2) for the pos.");
      return NULL;
    }
  
  //error check blockname to be the correct size
  if(strlen(blocknamein) > 4)
    {
      PyErr_Format(PyExc_TypeError,"'blockname' must be a four or less character string!");
      return NULL;
    }
  else
    {
      //put into blockname and make sure it is upper case
      for(i=0;i<strlen(blocknamein);++i)
	blockname[i] = blocknamein[i];
      for(i=strlen(blocknamein);i<4;++i)
	blockname[i] = ' ';
      for(i=0;i<4;++i)
	blockname[i] = toupper(blockname[i]);
    }
  
  //error check pos
  if(strcmp(blockname,"POS ") == 0 || strcmp(blockname,"VEL ") == 0 || strcmp(blockname,"ACCE") == 0)
    {
      if(pos < 0 || pos > 2)
	{
	  PyErr_Format(PyExc_TypeError,"pos has to be 0, 1, or 2 for reading positions, velocities, or accelerations!");
	  return NULL;
	}
    }
  
  //read from file(s) in a loop
  //***keep track of global particle index through pc and local file index through pc_new***
  for(i=0,pc=0;i<nfiles;i++,pc=pc_new)
    {
      if(nfiles > 1)  //multiple file names are different
	sprintf(buf,"%s.%d",file,i);
      else
	sprintf(buf,"%s",file);

      if(!(fd=fopen(buf,"r")))  //error check and complain
	{
	  PyErr_Format(PyExc_IOError,"gadgetPyIO module can't open file: '%s'",buf);
	  return NULL;
	}

      printf("reading: '%s' with blockname: '%s'\n",buf,blockname); 
      fflush(stdout);
      
      fread(&dummy,sizeof(dummy),1,fd);
      fread(&fileblockname,4*sizeof(char),1,fd);
      fseek(fd,dummy-4*sizeof(char),SEEK_CUR);
      fread(&dummy,sizeof(dummy),1,fd);
      
      if(strcmp(fileblockname,"HEAD") == 0)
	{
	  fread(&dummy, sizeof(dummy), 1, fd);
	  fread(&header, sizeof(header), 1, fd);
	  fread(&dummy, sizeof(dummy), 1, fd);
	}
      else
	{
	  PyErr_Format(PyExc_IOError,"first block is not labeled 'HEAD'!");
	  fclose(fd);
	  return NULL;
	}
      
      //allocate data array and get total number of particles to start reading
      if(i == 0)
	{
	  if(nfiles == 1)
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npart[k];
	      Ngas = header.npart[0];
	    }
	  else
	    {
	      for(k=0,NumPart=0;k<6;k++)
		NumPart += header.npartTotal[k];
	      Ngas = header.npartTotal[0];
	    }
	  
	  //make numpy array and copy info into it
	  dim[0] = NumPart;
	  if(strcmp(blockname,"ID  ") == 0)
	    array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_INT);
	  else
	    array = (PyArrayObject *) PyArray_SimpleNew(ndim, dim, NPY_FLOAT);
	}
      
      for(k=0,ntot_withmasses=0;k<6;k++)
	if(header.mass[k] == 0)
	  ntot_withmasses += header.npart[k];
      
      //loop through (and skip) blocks looking for the correct one
      fread(&dummy,sizeof(dummy),1,fd);
      fread(&fileblockname,4*sizeof(char),1,fd);
      fseek(fd,dummy-4*sizeof(char),SEEK_CUR);
      fread(&dummy,sizeof(dummy),1,fd);
      while(strcmp(fileblockname,blockname) != 0 && !feof(fd))
	{
	  fread(&dummy, sizeof(dummy), 1, fd);
	  fseek(fd,dummy,SEEK_CUR);
	  fread(&dummy, sizeof(dummy), 1, fd);
	  
	  fread(&dummy,sizeof(dummy),1,fd);
	  fread(&fileblockname,4*sizeof(char),1,fd);
	  fseek(fd,dummy-4*sizeof(char),SEEK_CUR);
	  fread(&dummy,sizeof(dummy),1,fd);
	}
      
      if(feof(fd) && strcmp(blockname,"MASS") == 0 && ntot_withmasses == 0)
	{
	  //assign masses from header
	  for(k=0,pc_new=pc;k<6;k++)
	    {
	      for(n=0;n<header.npart[k];n++)
		{
		  *((float*)PyArray_GETPTR1(array,pc_new)) = header.mass[k];
		  pc_new++;
		}
	    }
	}
      else if(feof(fd) && strcmp(blockname,"MASS") == 0 && ntot_withmasses > 0)
	{
	  //file structure is not consistent with a gadget2 snapshot - return and error
	  fclose(fd);
	  PyErr_Format(PyExc_IOError,"file structure information on particle masses is not consistent with a Gadget2 snapshot!");
	  Py_DECREF(array);
	  return NULL;
	}
      else if(feof(fd) && !(strcmp(blockname,"MASS") == 0))
	{
	  //complain a bit
	  fprintf(stdout,"could not find block '%s' in file '%s'!\n",blockname,buf);
	}
      else  //read the block!
	{
	  FLAG_NOBLOCK = 0;
	  fread(&dummy, sizeof(dummy), 1, fd);
	  
	  if(strcmp(blockname,"POS ") == 0 || strcmp(blockname,"VEL ") == 0 || strcmp(blockname,"ACCE") == 0)
	    {
	      for(k=0,pc_new=pc;k<6;k++)
		{
		  for(n=0;n<header.npart[k];n++)
		    {
		      switch(pos)
			{
			case 0:
			  fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
			  fseek(fd,2*sizeof(float),SEEK_CUR);
			  break;
			case 1:
			  fseek(fd,sizeof(float),SEEK_CUR);
			  fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
			  fseek(fd,sizeof(float),SEEK_CUR);
			  break;
			case 2:
			  fseek(fd,2*sizeof(float),SEEK_CUR);
			  fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
			  break;
			}
		      pc_new++;
		    }
		}
	    }
	  else if(strcmp(blockname,"ID  ") == 0)
	    {
	      for(k=0,pc_new=pc;k<6;k++)
		{
		  for(n=0;n<header.npart[k];n++)
		    {
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(int),1,fd);
		      pc_new++;
		    }
		}
	    }
	  else if(strcmp(blockname,"MASS") == 0)
	    {
	      for(k=0,pc_new=pc;k<6;k++)
		{
		  for(n=0;n<header.npart[k];n++)
		    {
		      if(header.mass[k] == 0)
			fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      else
			*((float*)PyArray_GETPTR1(array,pc_new)) = header.mass[k];
		      pc_new++;
		    }
		}
	    }
	  else if(strcmp(blockname,"POT ") == 0)
	    {
	      for(k=0,pc_new=pc;k<6;k++)
		{
		  for(n=0;n<header.npart[k];n++)
		    {
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      pc_new++;
		    }
		}
	    }
	  else if(strcmp(blockname,"TSTP") == 0)
	    {
	      for(k=0,pc_new=pc;k<6;k++)
		{
		  for(n=0;n<header.npart[k];n++)
		    {
		      fread(PyArray_GETPTR1(array,pc_new),sizeof(float),1,fd);
		      pc_new++;
		    }
		}
	    }
	  else if(strcmp(blockname,"U   ") == 0)
	    {
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  pc_sph++;
		}
	    }
	  else if(strcmp(blockname,"RHO ") == 0)
	    {
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  pc_sph++;
		}
	    }
	  else if(strcmp(blockname,"HSML") == 0)
	    {
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  pc_sph++;
		}
	    }
	  else if(strcmp(blockname,"ENDT") == 0)
	    {
	      for(n=0,pc_sph=pc;n<header.npart[0];n++)
		{
		  fread(PyArray_GETPTR1(array,pc_sph),sizeof(float),1,fd);
		  pc_sph++;
		}
	    }
	}
      
      //increase pc_new to keep indexing working
      for(k=0,pc_new=pc;k<6;k++)
	pc_new += header.npart[k];
      
      fclose(fd);
    }
  
  if(FLAG_NOBLOCK && !(strcmp(blockname,"MASS") == 0))
    {
      fprintf(stdout,"could not find block '%s' in any file!\n",blockname);
      Py_DECREF(array);
      Py_RETURN_NONE;
    }
  else
    {
      array = (PyArrayObject *) PyArray_Return(array);
      return (PyObject *) array;
    }
}

PyMethodDef methods[] = {
  {"readOne", readone, METH_VARARGS, "for all particles read a single phase space variable: x = gadgetPyIO.readOne(filename,nfiles,0), y = gadgetPyIO.readOne(filename,nfiles,1) \n vx = gadgetPyIO.readOne(filename,nfiles,0), vy = gadgetPyIO.readOne(filename,nfiles,1)"},
  {"readIds", readids, METH_VARARGS, "for all particles read their ids: ids = gadgetPyIO.readIds(filename,nfiles) "},
  {"readTypes", readparttype, METH_VARARGS, "for all particles read their types: type = gadgetPyIO.readTypes(filename,nfiles) "},
  {"readTypesTypeTwo", readparttypetypetwo, METH_VARARGS, "for all particles read their types from a type 2 snapshot: type = gadgetPyIO.readTypesTypeTwo(filename,nfiles) "},
  {"readMasses", readmasses, METH_VARARGS, "for all particles read their masses: m = gadgetPyIO.readMasses(filename,nfiles) "},
  {"printHeaders", printHeaders, METH_VARARGS, "for all files print the headers: gadgetPyIO.printHeaders(filename,nfiles) "},
  {"printHeadersTypeTwo", printHeadersTypeTwo, METH_VARARGS, "for all type 2 files print the headers: gadgetPyIO.printHeaders(filename,nfiles) "},
  {"readGasProps", readgasprops, METH_VARARGS, "for gas particles read their SPH properties: prop = gadgetPyIO.readGasProps(filename,nfiles,'prop',[DOCONV,UnitLength_in_cm,UnitMass_in_g,UnitVelocity_in_cm_per_s]) \n prop can be one of 'RHO ' 'U   ' or 'HSML' \n If DOCONV is true (i.e. not zero) its units are converted according the passed internal units or using the preferred internal units of Gadget2 (h^{-1}kpc, 1e10 h^{-1}(Solar masses), km/s) "},
  {"readBlockTypeTwo", readBlockTypeTwo, METH_VARARGS, "for a type 2 snapshot file, read a block of data specified by name via x = gadgetPyIO.readBlockTypeTwo(<filename>,nfile,'blockname',[pos]) \n If the block name is 'POS ', 'VEL ', or 'ACCE' pos must be an integer corresponding to (x,y,z) = (0,1,2)"},
    {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC 
initgadgetPyIO()
{
    (void) Py_InitModule("gadgetPyIO", methods);   
    import_array();
}

/* this template shows how one may convert from Gadget's units
 * to cgs units.
 * In this example, the temperate of the gas is computed.
 * (assuming that the electron density in units of the hydrogen density
 * was computed by the code. This is done if cooling is enabled.)
 int unit_conversion(void)
 {
 double GRAVITY, BOLTZMANN, PROTONMASS;
 double UnitLength_in_cm, UnitMass_in_g, UnitVelocity_in_cm_per_s;
 double UnitTime_in_s, UnitDensity_in_cgs, UnitPressure_in_cgs, UnitEnergy_in_cgs;  
 double G, Xh, HubbleParam;
 
 int i;
 double MeanWeight, u, gamma;
 
 //physical constants in cgs units
 GRAVITY   = 6.672e-8;
 BOLTZMANN = 1.3806e-16;
 PROTONMASS = 1.6726e-24;
 
 //internal unit system of the code
 UnitLength_in_cm= 3.085678e21;   //code length unit in cm/h
 UnitMass_in_g= 1.989e43;         //code mass unit in g/h
 UnitVelocity_in_cm_per_s= 1.0e5;
 
 UnitTime_in_s= UnitLength_in_cm / UnitVelocity_in_cm_per_s;
 UnitDensity_in_cgs= UnitMass_in_g/ pow(UnitLength_in_cm,3);
 UnitPressure_in_cgs= UnitMass_in_g/ UnitLength_in_cm/ pow(UnitTime_in_s,2);
 UnitEnergy_in_cgs= UnitMass_in_g * pow(UnitLength_in_cm,2) / pow(UnitTime_in_s,2);
 
 G=GRAVITY/ pow(UnitLength_in_cm,3) * UnitMass_in_g * pow(UnitTime_in_s,2);
 
 
 Xh= 0.76;  //mass fraction of hydrogen
 HubbleParam= 0.65;
 
 for(i=1; i<=NumPart; i++)
 {
 if(P[i].Type==0)  //gas particle
 {
 MeanWeight= 4.0/(3*Xh+1+4*Xh*P[i].Ne) * PROTONMASS;
 
 //convert internal energy to cgs units
 
 u  = P[i].U * UnitEnergy_in_cgs/ UnitMass_in_g;
 
 gamma= 5.0/3;
 
 //get temperature in Kelvin
 
 P[i].Temp= MeanWeight/BOLTZMANN * (gamma-1) * u;
 }
 }
 }
*/
